'''
Created on Mar 3, 2011

@author: oabalbin
'''
import os
import sys
import pickle
import numpy as np

from optparse import OptionParser
from collections import defaultdict, deque
from exome.plotting.plots import plot_scatter2, plot_hist, plot_hist2, plot_mat2d, plot_hist_ratio
from exome.vcftools.vcf import *


class commonSNVs:
    def __init__(self, vcf_keys, vcf_a1,vcf_a2, uniqeKeys1=[], unique2vcf1=[], uniqeKeys2=[], unique2vcf2=[]):
        self.snps=vcf_keys
        self.vcf1=vcf_a1
        self.vcf2=vcf_a2
        self.uvcf1=unique2vcf1
        self.uvcf2=unique2vcf2
        self.ukeys1=uniqeKeys1
        self.ukeys2=uniqeKeys2


def exctract_info(vcffile, tags, caller):
    '''
    Returns to numpy arrays one with keys, 
    one n-dim with the info on it 
    '''
    vcf=VCF()
    vcf.readFromFile(vcffile)
    #tags = ["QUAL","DP","_SAMPLES","ADref","ADalt"]
    
    # arrays
    vcf_keys=np.empty(len(vcf.keys()),dtype='str')
    vcf_a=np.zeros( (len(vcf.keys()),len(tags)) )
    
    for i, snp in enumerate(vcf):
        vcf_keys[i]=snp['COORD']
        vcf_a[i,0]=snp[tags[0]]
        vcf_a[i,1]=snp[tags[1]]
        # Returns a set with two dictionaries
        vcf_a[i,2],vcf_a[i,3], vcf_a[i,4]= get_fraction_variant(snp,tags, caller)
        
    tm= commonSNVs(vcf_keys[:i], vcf_a[:i,:],[], 
               [], [], [], [])
        
    return tm

def get_fraction_variant(snp,tags, caller):
    '''
    For Sam tools
    Number of 1) forward ref alleles; 
    2) reverse ref; 3) forward non-ref; 
    4) reverse non-ref alleles, used in variant calling. 
    Sum can be smaller than DP because low-quality bases are not counted.
    '''
    sam_cov="DP4"
    ZERO=-1.0
    ref, alt, frac=0.0,0.0,0.0
    samples = snp[tags[2]]
    if caller=="gatk":
        for sp in samples:
            ref += int(sp[tags[3]])
            alt += int(sp[tags[4]])    
        frac = float(alt)/(ref+alt)# if ref !=0.0 else ZERO
    
    elif caller == "sam":
        counts = snp[sam_cov].split(',')
        ref = int(counts[0]) + int(counts[1]) 
        alt= int(counts[2]) + int(counts[3])
        frac = float(alt)/(ref+alt)# if ref !=0.0 else ZERO
        
    
    return ref, alt, frac
    
        


def get_qual_dp(sites, vcfObj):
    '''
    '''
    tags = ["QUAL","DP","INDEL", "InDel"]
    indelct=0
    vcf_a1=np.zeros( (len(sites),3) )
    vcf_keys=np.empty(len(sites),dtype='str')
    i=0
    for j, pos in enumerate(sites):
        thsnps1= vcfObj.get([pos])
            
        for snp1 in thsnps1: 
            if (tags[2] in snp1.keys() or tags[3] in snp1.keys()): # In GATK the InDel tag is associated to the FILTER filed and can e like InDel;STAND_FILTER
                indelct += 1         
                continue
            
            vcf_keys[i]=pos
            vcf_a1[i,0],vcf_a1[i,1]=snp1[tags[0]], snp1[tags[1]]
            i+=1

    print >> sys.stderr,"The following number of indels was excluded from further analysis %d"%indelct
    
    return vcf_keys[:i], vcf_a1[:i,:]



def exctract_common_snps_info(vcffile1, vcffile2, tags, ftags):
    '''
    '''
    # Read vcf files
    vcf1=VCF()
    vcf1.readFromFile(vcffile1)
    vcf2=VCF()
    vcf2.readFromFile(vcffile2)
    # Get snps positions
    vcf_keys1 = vcf1.keys()
    vcf_keys2 = vcf2.keys()
    # Get common sites
    common_sites = set(vcf_keys1).intersection(set(vcf_keys2)) 
    
    unique_vcf1=set(vcf_keys1).difference(set(vcf_keys2))
    unique_vcf2=set(vcf_keys2).difference(set(vcf_keys1))
    ukeys1,uvcf_a1 = get_qual_dp(unique_vcf1, vcf1)
    ukeys2,uvcf_a2 = get_qual_dp(unique_vcf2, vcf2)
    
    ## Get Information
    #tags = ["QUAL","DP","INDEL", "InDel"]
    vcf_keys=np.empty(len(common_sites),dtype='str')
    vcf_a1=np.zeros( (len(common_sites),len(tags)) )
    vcf_a2=np.zeros( (len(common_sites),len(tags)) )
    indelct=0
    i=0
    for j, pos in enumerate(common_sites):
        thsnps1= vcf1.get([pos])
        thsnps2= vcf2.get([pos])
        if len(thsnps1) != len(thsnps2):
            print >> sys.stderr,"Problem with the snps dimensions between vcf1, vcf2"
        
        for snp1, snp2 in zip(thsnps1,thsnps2): 
            if (ftags[0] in snp1.keys() or ftags[0] in snp2.keys()) or \
               (ftags[1] in snp1.keys() or ftags[1] in snp2.keys()): # In GATK the InDel tag is associated to the FILTER filed and can e like InDel;STAND_FILTER
                indelct += 1
                continue
                
            vcf_keys[i]=pos
            vcf_a1[i,0],vcf_a1[i,1]=snp1[tags[0]], snp1[tags[1]]
            vcf_a2[i,0], vcf_a2[i,1]=snp2[tags[0]], snp2[tags[1]]
            i+=1
    
    # class with common and unique SNVs
    tm= commonSNVs(vcf_keys[:i], vcf_a1[:i,:],vcf_a2[:i,:], 
                   ukeys1, uvcf_a1, ukeys2, uvcf_a2)

    print >> sys.stderr,"The following number of indels was excluded from further analysis %d"%indelct
    
    return tm #vcf_keys, vcf_a1, vcf_a2

def get_SNVcall_prob(qual):
    '''
    Calculate the probability of SNV
    call based on the phred score
    '''
    return  1 - 10**((-1*qual)/10.0)

def get_SNVcall_zcore(qual):
    '''
    normalize the QUAL phred values
    '''
    return (qual - np.mean(qual))/np.std(qual)
    
def get_SNVcall_copa(qual):
    '''
    Normalize using the copa score of the QUAL
    '''
    mad = np.median(np.abs(qual - np.median(qual)))
        
    cqual = (qual - np.median(qual))/mad
    return cqual


def create_qual_isec_mat(qual1, qual2, dbins):
    '''
    All bins in the histogram equal width
    '''
    if len(qual1) != len(qual2):
        print >> sys.stderr,"Different dimmension for common snps between vcf1, vcf2"
        sys.exit(0)
    
    ind = np.arange(0,len(qual1))

    # make histogram for quala and qualb
    h1, db1 = np.histogram(qual1,dbins)
    h2, db2 = np.histogram(qual2,dbins)
    
    # Tells what bin each element of qual belongs to
    hind1, hind2 = np.digitize(qual1, db1), np.digitize(qual2, db2)
    
    # Find a set for each digitized values
    sind1, sind2 =set(hind1),set(hind2) 
    
    # loop between the values
    qual_mat=np.zeros((len(db1),len(db2)),dtype='int')
    
    for i in sind1:
        thelem_i = ind[hind1==i]
        for j in sind2:
            thelem_j = ind[hind2==j]
            ct = len(set(thelem_i).intersection(set(thelem_j)))
            # i-1,j-1 because in the binning the numbers start at 1
            # instead of 0
            qual_mat[i-1,j-1] = ct
    
    tm = type('qualmat', (object,),
         dict(mat=qual_mat,
              xedges=db1,
              yedges=db2))
    
    return tm
        

def create_qual_isec_mat2(qual1, qual2, dbins, pbins, phred, dmin,dmax):
    '''
    Histogram with bins of different width.
    It rescales different any values > phred
    according to pbins
    '''
    if len(qual1) != len(qual2):
        print >> sys.stderr,"Different dimmension for common snps between vcf1, vcf2"
        sys.exit(0)
    
    d = np.linspace(dmin,0.9,dbins)
    p = np.linspace(phred,1,pbins)
    dpbins = np.append(d,p)
    ind = np.arange(0,len(qual1))
    

    # make histogram for quala and qualb
    h1, db1 = np.histogram(qual1,dpbins)
    h2, db2 = np.histogram(qual2,dpbins)
    print db1
    print dpbins
    # Tells what bin each element of qual belongs to
    hind1, hind2 = np.digitize(qual1, db1), np.digitize(qual2, db2)
    
    # Find a set for each digitized values
    sind1, sind2 =set(hind1),set(hind2) 
    
    # loop between the values
    qual_mat=np.zeros((len(db1),len(db2)),dtype='int')
    
    for i in sind1:
        thelem_i = ind[hind1==i]
        for j in sind2:
            thelem_j = ind[hind2==j]
            ct = len(set(thelem_i).intersection(set(thelem_j)))
            # i-1,j-1 because in the binning the numbers start at 1
            # instead of 0
            qual_mat[i-1,j-1] = ct
    
    tm = type('qualmat', (object,),
         dict(mat=qual_mat,
              xedges=db1,
              yedges=db2))
    
    return tm


def create_qual_isec_mat3(qual1, qual2, uqual1, uqual2, dbins, pbins, phred, dmin,dmax):
    '''
    It plots shared and unique positions.
    uquala, uqualb quality score of positions not shared
    between arrays
    '''
    if len(qual1) != len(qual2):
        print >> sys.stderr,"Different dimmension for common snps between vcf1, vcf2"
        sys.exit(0)
    
    d = np.linspace(dmin,0.9,dbins)
    p = np.linspace(phred,1,pbins)
    dpbins = np.append(d,p)
    ind = np.arange(0,len(qual1))
    uind1 = np.arange(0,len(uqual1))
    uind2 = np.arange(0,len(uqual2))

    # make histogram for quala and qualb, uquala, uqualb
    h1, db1 = np.histogram(qual1,dpbins)
    h2, db2 = np.histogram(qual2,dpbins)
    h3, db3 = np.histogram(uqual1,dpbins)
    h4, db4 = np.histogram(uqual2,dpbins)

    
    # Tells what bin each element of qual belongs to
    hind1, hind2 = np.digitize(qual1, db1), np.digitize(qual2, db2)
    hind3, hind4 = np.digitize(uqual1, db3), np.digitize(uqual2, db4)    
    # Find a set for each digitized values
    sind1, sind2 =set(hind1),set(hind2) 
    # Find a set for each digitized values
    sind3, sind4 =set(hind3),set(hind4)
    
    # loop between the values
    qual_mat=np.zeros((len(db1),len(db2)),dtype='int')
    
    for i in sind1:
        thelem_i = ind[hind1==i]
        for j in sind2:
            thelem_j = ind[hind2==j]
            ct = len(set(thelem_i).intersection(set(thelem_j)))
            # i-1,j-1 because in the binning the numbers start at 1
            # instead of 0
            qual_mat[i-1,j-1] = ct
    
    # Include the SNVs called for one but 
    #not the other method.
    j=1
    for i in sind3:
        thelem_i = len(uind1[hind3==i])
        qual_mat[i-1,j-1] += thelem_i
    i=1
    for j in sind4:
        thelem_j = len(uind2[hind4==j])
        qual_mat[i-1,j-1] += thelem_j
        
    tm = type('qualmat', (object,),
         dict(mat=qual_mat,
              xedges=db1,
              yedges=db2))
    
    return tm



    

def read_files(vcffile1, vcffile2, fname, caller):
    '''
    '''
    # Do not change the order of these tags
    tags = ["QUAL","DP","_SAMPLES","ADref","ADalt"]
    ftags = ["INDEL", "InDel"]
    
    if vcffile2 is not None:
        tm =  exctract_common_snps_info(vcffile1, vcffile2, tags, ftags)
    else:
        tm =  exctract_info(vcffile1, tags, caller)

    of=open(fname,'w')
    pickle.dump(tm,of)
    return tm


def bin_size_1d(mat1, method):
    '''
    Scott, 1979
    W = 3.49*sd*N^(-1/3)
    Izemman 1991
    W = 2*(IQR)*N*(-1/3)
    '''
    from scipy import stats
    m75 = stats.scoreatpercentile(mat1, 75)
    m25 = stats.scoreatpercentile(mat1, 25)
    msd = np.std(mat1)
    N = len(mat1)
    if method=='scott':
        W = 2*(msd)*(N**(-1.0/3.0))
    elif method == 'izemman':
        W = 2*(m75 - m25)*(N**(-1.0/3.0))
    
    print W, N, m75, m25, msd
    return W

def bins_edges(mat1,bwidth):
    '''
    Create the bin edges
    '''
    mmin, mmax = min(mat1), max(mat1)
    ebins = np.arange(mmin,mmax,bwidth)
    return ebins
    

def main_compare(comSNVs,folderpath, dbSNPsSNVs=None):
    '''
    '''
    # dictionary of fields in comSNVs mat
    pd ={'qual':0,'DP':1}
    DPlabels={'xlabel':'Depth of coverage','ylabel':'SNV call Quality score'}
    DPlabels2={'xlabel':'Depth of coverage GATK','ylabel':'Depth of coverage Samtools'}
    DPlabels3={'xlabel':'Quality score GATK','ylabel':'Quality score Samtools'}
    DPlabels4={'xlabel':'Quality score','ylabel':'Ratio HapMap (TP)/Total SNV calls'}
    H2labels={'xlabel':'SamTools SNV Quality score','ylabel':'GATK SNV Quality score'}
    
    
    
    histname2d='hist2d_Qscore'
    histname3d='hist3d_Qscore'
    histname1='hist_GATK_Qscore_common'
    histname2='hist_SAM_Qscore_common'
    histname1b='hist_GATK_Qscore_all'
    histname2b='hist_SAM_Qscore_all'
    histname1c='hist_GATK_HDscore_all'
    histname2c='hist_SAM_HDscore_all'
    histratio1b='hist_ratio_GATK_TP'
    histratio2b='hist_ratio_SAM_TP'
    
    ####
    histname3='hist_GATK_QscoreZ_common'
    histname4='hist_SAM_QscoreZ_common'
    histname2df='hist2d_Qscore_full'
    histname3df='hist3d_Qscore_full'
    
    scatter1='GATK_qual_DP'
    scatter2='SAM_qual_DP'
    scatter3='GATK_SAM_DP'
    scatter4='GATK_SAM_Qscore'
    
    matf='hist_QSCORE.txt'
 
    dbins=10
    pbins=10
    hbins=100
    hdbins=50
    phred=0.9999 # 40
    dmin,dmax = 0.0, 1.0
    ##########
    print len(comSNVs.snps)
    ##########
    qual1 = get_SNVcall_prob(comSNVs.vcf1[:,pd['qual']])
    qual2 = get_SNVcall_prob(comSNVs.vcf2[:,pd['qual']])
    
    ######## Plot Matrices
    mat2d = create_qual_isec_mat(qual1, qual2,dbins)
    plot_mat2d(mat2d.mat,mat2d.xedges, mat2d.yedges, H2labels, histname2d, folderpath)
    #print mat2d.mat
    mat3d = create_qual_isec_mat2(qual1, qual2,dbins,pbins, phred, dmin,dmax)
    plot_mat2d(mat3d.mat,mat3d.xedges, mat3d.yedges, H2labels, histname3d, folderpath)
    
    ### Test (qual1, qual2, uqual1, uqual2, dbins, pbins, phred, dmin,dmax)
    uqual1 = get_SNVcall_prob(comSNVs.uvcf1[:,pd['qual']])
    uqual2 = get_SNVcall_prob(comSNVs.uvcf2[:,pd['qual']])
    mat4d = create_qual_isec_mat3(qual1, qual2, uqual1, uqual2, dbins,pbins, phred, dmin,dmax)
    plot_mat2d(mat4d.mat,mat4d.xedges, mat4d.yedges, H2labels, histname3df, folderpath)

    ouf=open(os.path.join(folderpath,matf),'w')
    for i in range(mat3d.mat.shape[0]):
        ouf.write(",".join(map(str,list(mat3d.mat[i,:]))).replace(',','\t')+'\n')
    
    #########
    ######### Plot histograms
    
    #plot_hist(vcf_a1_qprob,histname1,folderpath)
    #plot_hist(vcf_a2_qprob,histname2,folderpath)
    #### Plot common snps
    plot_hist(comSNVs.vcf1[:,pd['qual']],histname1,folderpath,hbins)
    plot_hist(comSNVs.vcf2[:,pd['qual']],histname2,folderpath,hbins)    
    hqual3=get_SNVcall_copa(comSNVs.vcf1[:,pd['qual']])
    hqual4=get_SNVcall_copa(comSNVs.vcf2[:,pd['qual']])
    plot_hist(hqual3,histname3,folderpath,hbins)
    plot_hist(hqual4,histname4,folderpath,hbins)
    
    #### Plot all SNVs common and not shared
    if dbSNPsSNVs is not None:
        allgatk = np.append( comSNVs.vcf1[:,pd['qual']], comSNVs.uvcf1[:,pd['qual']])
        allsam = np.append( comSNVs.vcf2[:,pd['qual']], comSNVs.uvcf2[:,pd['qual']])
        dbsnpgatk=np.append( dbSNPsSNVs.vcf1[:,pd['qual']], dbSNPsSNVs.uvcf1[:,pd['qual']])
        dbsnpsam=np.append( dbSNPsSNVs.vcf2[:,pd['qual']], dbSNPsSNVs.uvcf2[:,pd['qual']])
        
        hdbins_gatk = bins_edges(allgatk, bin_size_1d(allgatk, 'scott'))
        hdbins_sam = bins_edges(allsam, bin_size_1d(allsam, 'scott'))
        print len(hdbins_gatk), len(hdbins_sam)
        
        plot_hist2(allgatk, dbsnpgatk, histname1b, folderpath,hbins)#hdbins_gatk)
        plot_hist2(allsam, dbsnpsam, histname2b, folderpath, hdbins_sam)
        
        plot_hist_ratio(allgatk, dbsnpgatk, histratio1b, folderpath,hdbins_gatk,DPlabels4) # bins_edges(allgatk,hbins)
        plot_hist_ratio(allsam, dbsnpsam, histratio2b, folderpath, hdbins_sam, DPlabels4)
    
    
    
    #### Plot QD for all SNVs common and not shared
    if dbSNPsSNVs is not None:        
        allgatk = np.true_divide(np.append( comSNVs.vcf1[:,pd['qual']], comSNVs.uvcf1[:,pd['qual']]),
                                 np.append( comSNVs.vcf1[:,pd['DP']], comSNVs.uvcf1[:,pd['DP']]))
        allsam = np.true_divide(np.append( comSNVs.vcf2[:,pd['qual']], comSNVs.uvcf2[:,pd['qual']]),
                                np.append( comSNVs.vcf2[:,pd['DP']], comSNVs.uvcf2[:,pd['DP']]))
        dbsnpgatk = np.true_divide(np.append( dbSNPsSNVs.vcf1[:,pd['qual']], dbSNPsSNVs.uvcf1[:,pd['qual']]),
                                   np.append( dbSNPsSNVs.vcf1[:,pd['DP']], dbSNPsSNVs.uvcf1[:,pd['DP']]))
        dbsnpsam= np.true_divide(np.append( dbSNPsSNVs.vcf2[:,pd['qual']], dbSNPsSNVs.uvcf2[:,pd['qual']]),
                                 np.append( dbSNPsSNVs.vcf2[:,pd['DP']], dbSNPsSNVs.uvcf2[:,pd['DP']]))
        
        hdbins_gatk = bins_edges(allgatk, bin_size_1d(allgatk, 'scott'))
        hdbins_sam = bins_edges(allsam, bin_size_1d(allsam, 'scott'))
        print len(hdbins_gatk), len(hdbins_sam)
        
        plot_hist2(allgatk, dbsnpgatk, histname1c, folderpath, hdbins_gatk)
        plot_hist2(allsam, dbsnpsam, histname2c, folderpath, hdbins_sam)

    
    
    ########
    plot_scatter2(comSNVs.vcf1[:,pd['DP']], comSNVs.vcf1[:,pd['qual']],
                  scatter1,folderpath,DPlabels)
    plot_scatter2(comSNVs.vcf2[:,pd['DP']], comSNVs.vcf2[:,pd['qual']],
                  scatter2,folderpath,DPlabels)
    
    plot_scatter2(comSNVs.vcf1[:,pd['DP']], comSNVs.vcf2[:,pd['DP']],
                  scatter3,folderpath, DPlabels2)
    
    plot_scatter2(comSNVs.vcf1[:,pd['qual']], comSNVs.vcf2[:,pd['qual']],
                  scatter4,folderpath, DPlabels3)


    sys.exit(0)



def main_single_analysis(comSNVs,folderpath, dbSNPsSNVs=None):
    '''
    GT:AD:DP:GL:GQ  0/1:0,1:1:-4.26,-0.30,-0.00:1.76        0/0:5,0:4:-0.00,-1.21,-13.18:12.03
    '''
    
    pd ={'qual':0,'DP':1,'RefDP':2,'AltDP':3,'FRAG':4}
    DPlabels={'xlabel':'Depth of coverage','ylabel':'SNV call Quality score'}
    scatter1='low_qualHapMap_GATK_qual_DP'
    scatter2='zoom_low_qualHapMap_GATK_qual_DP'
    DPlabels3={'xlabel':'ALT depth of coverage','ylabel':'SNV call Quality score'}
    DPlabels4={'xlabel':'ALT Fraction of coverage','ylabel':'SNV call Quality score'}
    scatter3='low_qualHapMap_GATK_qual_AltDP'
    scatter4='low_qualHapMap_GATK_qual_FRAG' 
    # Plot depth of coverage, and QUAL score
    print len(comSNVs.vcf1[:,pd['DP']])
    print np.mean(comSNVs.vcf1[:,pd['DP']]), min(comSNVs.vcf1[:,pd['DP']]), max(comSNVs.vcf1[:,pd['DP']])
    print np.mean(comSNVs.vcf1[:,pd['qual']]),min(comSNVs.vcf1[:,pd['qual']]), max(comSNVs.vcf1[:,pd['qual']])
    plot_scatter2(comSNVs.vcf1[:,pd['DP']], comSNVs.vcf1[:,pd['qual']],
                  scatter1,folderpath,DPlabels)
    
    ind=np.arange(len(comSNVs.vcf1[:,pd['DP']]))
    dep_low = ind[comSNVs.vcf1[:,pd['DP']] < 100]
    qual_low = ind[comSNVs.vcf1[:,pd['qual']] < 100]
    ove_low = list(set(dep_low).intersection(set(qual_low)))
    plot_scatter2(comSNVs.vcf1[ove_low,pd['DP']], comSNVs.vcf1[ove_low,pd['qual']],
                  scatter2,folderpath,DPlabels)
    
    plot_scatter2(comSNVs.vcf1[:,pd['AltDP']], comSNVs.vcf1[:,pd['qual']],
                  scatter3,folderpath,DPlabels3)
    
    plot_scatter2(comSNVs.vcf1[:,pd['FRAG']], comSNVs.vcf1[:,pd['qual']],
                  scatter4,folderpath,DPlabels4)


    
    

    sys.exit(0)
    

if __name__ == '__main__':
        
    optionparser = OptionParser("usage: %prog [options] ")
    
    # Specific options
    optionparser.add_option("-q", "--query_vcf_file", dest="query_vcf_file",
                            help="query_vcf_file")
    optionparser.add_option("-o", "--other_vcf_file", dest="other_vcf_file",
                            help="others_vcf_file")
    optionparser.add_option("-f", "--outfolder", dest="outfolder",
                            help="outfolder")
    optionparser.add_option("-d", "--use_dump_files", dest="use_dump_files", 
                             action="store_true", default=False,
                            help="use dump files instead of reading the vfcfiles")
    optionparser.add_option("-b", "--use_dbSNP", dest="use_dbSNP", 
                             action="store_true", default=False,
                            help="use dbSNPdump files")
    optionparser.add_option("-p", "--use_dump_dbSNPfiles", dest="use_dump_dbSNPfiles", 
                            help="use use_dump_dbSNPfiles files instead of reading the vfcfiles")
    optionparser.add_option("-c", "--caller", dest="caller", 
                            help="caller used, gatk or sam")


    
    (options, args) = optionparser.parse_args()
    
    print options
    
    
    logfile = open(os.path.join(options.outfolder,'README.log'),'w')
    logfile.write('Run call:\n')
    [logfile.write(",".join([str(opt),str(val)]).replace(',','\t')+'\n') for opt, val in options.__dict__.iteritems()]
    
    if options.use_dump_files:
        commSNVs = pickle.load(open(os.path.join(options.outfolder,'commSNVs.dump')))
        
    if (options.other_vcf_file is not None) and (options.query_vcf_file):
        commSNVs = read_files(options.query_vcf_file, options.other_vcf_file,
                          os.path.join(options.outfolder,'commSNVs.dump'))
        
        if options.use_dbSNP:
            dbSNPsSNVs = pickle.load(open(options.use_dump_dbSNPfiles))
            main_compare(commSNVs,options.outfolder, dbSNPsSNVs)
        else:
            main_compare(commSNVs,options.outfolder)

    else:        
        commSNVs = read_files(options.query_vcf_file, options.other_vcf_file,
                          os.path.join(options.outfolder,'commSNVs.dump'),options.caller)
        if options.use_dbSNP:
            main_single_analysis(commSNVs,options.outfolder, dbSNPsSNVs)
        else:
            main_single_analysis(commSNVs,options.outfolder)
            
    
    
    

    

